TOI-2525 b and c: A Pair of Massive Warm Giant Planets with Strong Transit Timing Variations Revealed by TESS^*

Object TOI-2525 was observed in sectors 1, 3–11, and 13 during the first year of the TESS primary mission with a 30 minute cadence and sectors 27, 28, and 30–39 with a 2 minute cadence in the third mission year. Objects TOI-2525 b and c were identified in the light curves extracted from the TESS full-frame images (FFIs) using the tesseract ⁴¹  pipeline (F. Rojas et al. 2023, in preparation). A brief introduction of our FFI extraction with tesseract in the context of the WINE collaboration can be found in Schlecker et al. (2020), Gill et al. (2020), and Trifonov et al. (2021). Figure 1 shows the target pixel file image of TOI-2525 constructed from the TESS FFI image frames and Gaia DR3 data (Gaia Collaboration et al. 2021). Figure 1 shows that there are no bright contaminators in the FFI aperture (red solid contour); thus, we concluded that the transit signals are indeed coming from TOI-2525 and not from neighboring stars. This was later confirmed by ground-based transit detections of TOI-2525 b and c. However, a fainter source labeled "1" is occasionally in the FFI aperture, which dilutes the FFI light curves. Since tesseract does not correct contamination in the TESS apertures from nearby stars, we follow the same methodology as used in Trifonov et al. (2021) to calculate and apply a dilution correction for the contamination of TOI-2525 on the FFI light curves.

Zoom In Zoom Out Reset image size

We retrieved the 2 minute cadence light curves from the Mikulski Archive for Space Telescopes. ⁴² The Science Processing Operations Center (Jenkins et al. 2016) provides simple aperture photometry (SAP) and systematics-corrected presearch data conditioning (PDC) photometry (Smith et al. 2012; Stumpe et al. 2012). The PDCSAP light curves are corrected for contamination from nearby stars and instrumental systematics originating from, e.g., pointing drifts, focus changes, and thermal transients. In our work on TOI-2525, for the 2 minute cadence data, we only use the corrected PDCSAP data.

The Antarctica Search for Transiting ExoPlanets (ASTEP; Guillot et al. 2015) instrument is a 40 cm Newton telescope installed in 2010 at the Concordia station, located at −75.06°S, 123.3°E and an altitude of ∼3230 m. ASTEP is a robotic telescope dedicated to photometric observations of the fields of stars and their exoplanets.

Due to the extremely low data transmission rate at the Concordia station, the data are processed automatically on site using an IDL-based aperture photometry pipeline (Mékarnia et al. 2016). The raw light curves of up to 1000 stars of the field are transferred to Europe on a server in Roma, Italy, and then available for deeper analysis. These data files contain each star's flux computed through 10 fixed circular aperture radii, so that optimal calibrated light curves can be extracted.

Thanks to the accurate TTV model prediction constructed on the TESS data, we scheduled successful observations with ASTEP of TOI-2525 b and c. For TOI-2525 b, we detected a full and a partial transit event on the nights of UT 2021-09-17 and UT 2022-06-23, respectively, whereas for TOI-2525c, we observed two full transit events on the nights of UT 2021-04-15 and UT 2022-07-02 and two partial transit events on the nights of UT 2021-06-3 and UT 2021-09-10.

A partial transit of TOI-2525c was observed with the Siding Spring Observatory station of the Observatoire Moana telescope network (OM-SSO). The OM-SSO is an RC Optical Systems RCOS20 f8.1 telescope with a focal length of 3980 mm. It is equipped with an FLI Microline 16803 camera with 4k × 4k pixels of 9 with a pixel scale of 047 and a field of view of $30^{\prime} \times 30^{\prime}$. Observations were taken using an Astrodon Exoplanet (clear blue blocking) filter. Observations of TOI-2525c were performed on 2021 October 29 covering an ingress. The adopted exposure time was 147 s, and the airmass ranged from 1.15 to 2. The OM-SSO data were processed with a dedicated automated pipeline adapted from a version that was initially developed for obtaining differential photometry of Las Cumbres Observatory Global Telescope (LCOGT) light curves (Espinoza et al. 2019).

We observed one partial transit of TOI-2525b and two partial transits of TOI-2525c using the LCOGT 1.0 m network (Brown et al. 2013) nodes at the Cerro Tololo Inter-American Observatory (CTIO) and South Africa Astronomical Observatory (SAAO). The 1 m telescopes are equipped with 4096 × 4096 SINISTRO cameras having a pixel scale of 0389 pixel⁻¹, resulting in a field of view of 26' × 26'. The TOI-2525b transit was observed from SAAO on 2022 January 11 in the Sloan $g^{\prime}$ and $i^{\prime}$ filters using a 43 target aperture. Transits of TOI-2525c were observed twice from CTIO on 2021 December 18 and 2022 February 5 in the Sloan $g^{\prime}$ and $i^{\prime}$ filters, respectively, using 40–47 target apertures. The LCOGT data reduction and photometric measurements were performed using the AstroImageJ (Collins et al. 2017) software package.

Object TOI-2525 was monitored with the Planet Finder Spectrograph (PFS; Crane et al. 2006, 2008, 2010) installed at the 6.5 m Magellan/Clay telescope at Las Campanas Observatory. It was observed with the iodine gas absorption cell of the instrument at four different observing runs between 2020 November 3 and 2021 November 16, adopting an exposure time of 1200 s and using a 3 × 3 CCD binning mode to minimize read noise. Object TOI-2525 was also observed without the iodine cell in order to generate the template for computing the RVs, which were derived following the methodology of Butler et al. (1996). The mean uncertainty of the PFS RVs of TOI-2525 is 5.7 m s⁻¹. The PFS RVs are presented in Table 1.

2.5.1. FEROS

We identify the transit events of TOI-2525 b and c in the first-year TESS FFI light curves. Our preliminary transit characterization methodology includes detrending of the tesseract FFI light curves with a robust (iterative) Matérn Gaussian process (GP) kernel via the wotan package (see Hippke et al. 2019) and a transit signal search with the transitleastsquares (TLS; Hippke & Heller 2019) algorithm.

Figure 2 shows our TLS results on the combined TESS FFI light-curve data of TOI-2525. We first detect the stronger transit signal of TOI-2525c with a period of ≈49.3 days. We filter this signal by applying a Keplerian transit model with the obtained parameters from the TLS, and we seek additional transit signals in the model residuals. We clearly find the shallower but more frequent transit signal of the inner planet TOI-2525b at a period of P_b = 23.3 days. Figure 2 shows many other significant TLS peaks with a signal detection efficiency (see Hippke & Heller 2019) of >8.3 that are nothing more than the integer subharmonics of the transit signals.

Zoom In Zoom Out Reset image size

From the FFI data, we found that the light-curve mid-transit time of TOI-2525b and TOI-2525c are not linear in time and show strong deviations in the expected times of transit, i.e., TTVs.

The available Doppler data of TOI-2525 are too few for an independent RV validation of the two-planet system (see Table 1). Since we identified the transit events in 2019, we have made many attempts to collect precise spectroscopic data from the southern hemisphere. However, TOI-2525 is faint and requires significant observational efforts and excellent sky conditions. This, in combination with the COVID-19 pandemic closures of the ESO and Las Campanas observatories, prevented us from obtaining sufficient RV data. Nonetheless, we conclude that at this point, the PFS data alone are adequate for validation of the system when combined with the transit periods from TESS and can contribute to the planetary mass estimates when combined with the strong TTVs.

4.1.1. Extraction of TTVs

The TTVs were estimated by fitting all of the available photometric data using the exoplanet software package (Foreman-Mackey et al. 2021). We used the descriptive model TTVOrbit, which assumes Keplerian orbits for each planet but allows for the central time of each transit to be a free parameter in the model. The parameters of the model are the values of R_p for each planet, their impact parameters b, the transit times, the stellar mass M_⋆ and density ρ_⋆, quadratic limb-darkening coefficients u₁ and u₂ for each instrument used, and parameters to describe trends and correlations in the data. Regarding the latter, we adopt a linear model in time for the ground-based transits, which are often only partially observed, and a GP for each of the short- and long-cadence TESS data sets. The GP kernel adopted is a damped simple harmonic oscillator (Foreman-Mackey et al. 2017a) with Q = 1/3, variance σ_GP, and correlation length parameter ρ_GP.

Table 3 lists the resulting light-curve parameters, whereas Figure 3 shows the resulting light-curve models applied to all photometric data sets. Both TOI-2525 transit signals exhibit strong TTV libration. The values of the fitted transit times for each transit are listed in Table A2.

Zoom In Zoom Out Reset image size

Table 3. Planetary Radii, Orbital Inclinations, and Stellar Density and Mass Estimates of the TOI-2525 System Derived during TTV Extraction

Parameter	Median	σ
rb (RJup)	0.88	0.02
rc (RJup)	0.98	0.02
ib (deg)	89.31	0.03
ic (deg)	89.96	0.03
ρ⋆ (g cm−3)	2.14	0.04
M⋆ (M⊙)	0.85	0.05

Note. The remaining transit light-curve data parameter estimates are listed in Table A1, whereas the individual transit times for TOI-2525 b and c are listed in Table A2.

Download table as:  ASCIITypeset image

4.2.1. Joint RV and TTV Analysis

For the joint RV and TTV orbital analysis of the TOI-2525 system, we followed a similar route to the one we used in Trifonov et al. (2021) for the modeling of the TOI-2022 system. We refer the reader to that paper for a more detailed description of the chosen methodology. Briefly, for TOI-2525, we performed orbital fitting with the Exo-Striker exoplanet toolbox ⁴³  (Trifonov 2019) by adopting the self-consistent dynamical model on the extracted TTVs from Section 4.1.1. The TTV model can provide a relatively inexpensive orbital and dynamical solution to the system in terms of CPU time. However, it also comes with a severe ambiguity in eccentricity versus dynamical planetary mass (e.g., Lithwick et al. 2012; Dawson et al. 2021; Trifonov et al. 2021). Therefore, we fit the TTVs jointly with the RV data from PFS to better constrain the planetary dynamical masses and eccentricities. The RV model is intrinsic to Exo-Striker, whereas the TTV model is wrapped around the TTVfast package (Deck et al. 2014). The fitted parameters for each planet were the RV semiamplitude K (which is automatically converted to the dynamical planetary mass m_p ), orbital period P, eccentricity e, argument of periastron ω, and mean anomaly M₀. Our TTV+RV modeling scheme allows a difference in the planetary inclinations (Δi ≥0°). Thus, we also model the orbital inclinations i and the difference between the line of node ΔΩ. All of these parameters are valid for BJD = 2,458,333.52, which is an arbitrary epoch chosen slightly before the first transit event of the inner planet. The RV data offset and jitter ⁴⁴  parameters of PFS added two more free parameters.

Finally, our dynamical model, particularly the planetary masses, also depends on the stellar mass estimate of TOI-2525, for which we adopt a fixed value of 0.849 M_⊙. We set the time step in the dynamical model to dt = 0.02 days to assure sufficient orbital resolution and accuracy.

We ran a nested sampling (NS) scheme (Skilling 2004), which allowed us to efficiently explore the complex parameter space of osculating orbital elements and study the parameter posteriors and overall dynamics. Our NS run was performed with the dynesty sampler (Speagle 2020), which is integrated into Exo-Striker. We ran 100 "live points" per fitted parameter using the "Dynamic" NS scheme, focused on 100% posterior convergence instead of log evidence (see Speagle 2020 for details). Parameter priors were estimated by running several experimental NS runs and adopting a wide range of uniform parameter priors. After a few consecutive NS runs, we narrowed the adequate parameter space to be explored. We note that we adopted very narrow priors on the orbital inclinations to ensure that the TTV+RV model parameter space is consistent with the inclination estimates extracted from the TTVs. Thus, we eliminate configurations that could explain the RV data but would not lead to transit events (i.e., impact parameters b_b and b_c, which are inconsistent with the light-curve signal). The final adopted parameter priors are listed in Table 4.

Table 4. NS Priors, Posteriors, and the Optimum $-\mathrm{ln}{ \mathcal L }$ Orbital Parameters of the Two-planet System TOI-2525 Derived by Joint Dynamical Modeling of TTVs (TESS, ASTEP, SSO, and LCOGT) and RVs (FEROS and PFS)

	Median and 1σ		Max. $-\mathrm{ln}{ \mathcal L }$		Adopted Priors
Parameter	Planet b	Planet c	Planet b	Planet c	Planet b	Planet c
K (m s−1)	${7.1}_{-0.3}^{+0.3}$	${44.4}_{-1.5}^{+1.5}$	7.1	44.4	${ \mathcal U }$(5.0,30.00)	${ \mathcal U }$(20.0,60.0)
P (days)	${23.288}_{-0.002}^{+0.001}$	${49.260}_{-0.001}^{+0.001}$	23.289	49.260	${ \mathcal U }$(23.2,23.4)	${ \mathcal U }$(49.1,49.4)
e	${0.159}_{-0.007}^{+0.012}$	${0.152}_{-0.005}^{+0.006}$	0.171	0.160	${ \mathcal U }$(0.0,0.4)	${ \mathcal U }$(0.0,0.4)
ω (deg)	${346.3}_{-0.7}^{+0.7}$	${21.8}_{-0.8}^{+0.8}$	346.5	22.0	${ \mathcal U }$(0.0,360.0)	${ \mathcal U }$(0.0,360.0)
M0 (deg)	${120.8}_{-0.5}^{+0.6}$	${71.2}_{-1.0}^{+1.2}$	120.6	71.6	${ \mathcal U }$(0.0,360.00)	${ \mathcal U }$(0.0,360.00)
λ (deg)	${107.0}_{-0.8}^{+0.8}$	${93.4}_{-1.0}^{+0.9}$	107.1	93.1	(derived)	(derived)
i (deg)	${89.96}_{-0.07}^{+0.08}$	${89.99}_{-0.06}^{+0.08}$	89.97	89.99	${ \mathcal N }$(90.0,0.1)	${ \mathcal N }$(90.0,0.1)
Ω (deg)	0.0	${2.0}_{-1.2}^{+1.9}$	0.0	1.9	(fixed)	${ \mathcal N }$(0.0,15.0)
Δi (deg)	${2.0}_{-1.2}^{+1.9}$	...	1.9	...	(derived)	...
a (au)	${0.1511}_{-0.0025}^{+0.0025}$	${0.2491}_{-0.0042}^{+0.0041}$	0.1511	0.2491	(derived)	(derived)
m (MJup)	${0.088}_{-0.004}^{+0.005}$	${0.709}_{-0.034}^{+0.034}$	0.089	0.710	(derived)	(derived)
ρ (g cm−3)	${0.174}_{-0.015}^{+0.016}$	${1.014}_{-0.076}^{+0.084}$	0.174	1.014	(derived)	(derived)
RVoff. PFS (m s−1)	−28.5${}_{-7.4}^{+7.1}$	...	−32.2	...	${ \mathcal U }$(−300.00,100.0)	...
RVjit. PFS (m s−1)	${26.8}_{-5.4}^{+6.8}$	...	21.5	...	${ \mathcal J }$(0.0,50.0)	...

Note. The orbital elements are in the Jacobi frame and valid for epoch BJD = 2,458,333.52. The adopted priors are listed in the rightmost columns, and their meanings are ${ \mathcal U }$, uniform; ${ \mathcal N }$, Gaussian; and ${ \mathcal J }$, Jeffrey's (log-uniform). The derived planetary posterior parameters of a and m are calculated taking into account the stellar parameter uncertainties (see note in Table 2).

Download table as:  ASCIITypeset image

Figures 4 and 5 show the TTV and RV data together with the best-fit joint dynamical model of TOI-2525. The left and right panels of Figure 4 show the TTVs of TOI-2525b and c, respectively, fitted with the best-fit TTV model. Figure 5 shows the RV component of the best-fit model applied to the PFS RVs. Figure 5 also shows the two FEROS RVs fitted independently to the best fit with an optimized offset. The middle and right panels of Figure 5 show a phase-folded representation of the RV signals of TOI-2525b and c, respectively, modeled with an osculating period. The signal of the inner planet TOI-2525b is strongly overshadowed by the dominating signal of the outer one due to the large difference in K amplitudes (see Table 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The final posterior probability distributions with an RV linear trend are shown in Figure A1. Our final estimates for TOI-2525b and c lead to planetary orbital periods of P_b = ${23.288}_{-0.002}^{+0.001}$ and P_c = ${49.260}_{-0.001}^{+0.001}$ days, eccentricities of e_b = ${0.159}_{-0.007}^{+0.012}$ and e_c = ${0.152}_{-0.005}^{+0.006}$, and dynamical masses of m_b = ${0.088}_{-0.004}^{+0.005}$ and m_c = ${0.709}_{-0.034}^{+0.034}$ M_Jup. The mutual inclination is constrained to Δi = ${2\buildrel{\circ}\over{.} 0}_{-1.2}^{+1.9}$. Given the planetary radii obtained during the TTV extraction in Section 4.1.1, we derive a remarkably low density of ρ_b = 0.174${}_{-0.015}^{+0.016}$ g cm⁻³ for the inner planet and ρ_c = 1.014${}_{-0.076}^{+0.084}$ g cm⁻³ for the outer one. The full list of posterior and maximum $-\mathrm{ln}{ \mathcal L }$ (i.e., best-fit) estimates was derived from the joint TTV+RV model and listed in Table 4.

4.2.2. Joint RV and Photodynamical Analysis

Following Trifonov et al. (2021), we inspected the Hill (see Gladman 1993) and angular momentum deficit (AMD; see Laskar & Petit 2017) stability criteria of the TOI-2525 system. In terms of the classical Hill stability criterion, the TOI-2525 planetary system is predicted to be stable. We estimate a mutual Hill distance of ∼7.4 R_Hill,m, which is above the widely accepted distance of ∼3.5 R_Hill,m for the system to be considered Hill-stable. However, accounting for the estimated orbital eccentricities, semimajor axes, mutual orbital inclinations, and planetary masses of TOI-2525 b and c, the AMD criterion suggests that the TOI-2525 planetary system is unstable. The AMD criterion is very sensitive to eccentricities. Thus, the moderately eccentric orbits of the pair are the reason for the negative AMD result.

As discussed in Trifonov et al. (2021), the AMD and Hill stability criteria are only a proxy for long-term stability and do not account for the system's dynamics near MMRs. Therefore, to test the long-term dynamics and possible MMR in the system, we adopt exactly the same N-body numerical setup used in our recent analysis of TOI-2202 (Trifonov et al. 2021). This is adequate because TOI-2202 and TOI-2525 share somewhat similar physical and orbital characteristics. We refer the reader to Section 5 of Trifonov et al. (2021) for more details about the stability test performed here for TOI-2525. Briefly, we performed numerical integrations using a custom version of the Wisdom–Holman N-body algorithm (Wisdom & Holman 1991), which directly adopts and integrates the Jacobi orbital elements from the posterior orbital analysis. We tested the stability of the TOI-2525 system up to 1 Myr with a small time step of 0.02 days for 10,000 randomly chosen samples from the achieved orbital parameter posteriors from the TTV+RV dynamical modeling scheme. We automatically monitored the evolution of the planetary semimajor axes, eccentricities, secular apsidal angle Δω = ω_b − ω_c, and first-order 2:1 MMR angles θ₁ = λ_b − 2λ_c + ω_b and θ₂ = λ_b − 2λ_c + ω_c (where λ_b,c = M_b,c + ω_b,c is the mean longitude of planets b and c, respectively; see Lee 2004). These angles are important for libration in secular apsidal alignment or MMR.

We found that all 10,000 samples examined are stable for 1 Myr. Figure 6 presents the resulting posterior probability distribution of some of the important dynamical properties of the system, such as mean period ratio P_rat, mean eccentricities ${\hat{e}}_{{\rm{b}}}$ and ${\hat{e}}_{{\rm{c}}}$, peak-to-peak amplitudes ampl. e_b and ampl. e_c, dynamical masses of the planets m_b and m_c, and orbital semimajor axes a_b and a_c. Figure 6 shows that the period ratio evolution is bimodal, oscillating around either ∼2.127 or, more plausibly, ∼2.126, which is too far from the exact 2:1 period ratio for the system to librate in the 2:1 MMR. This is confirmed by the circulation of θ₁ and θ₂ (for which distributions are not shown in Figure 6). On the other hand, we observe a bimodality of the evolution of Δω, which is intriguing. About half of the sampled initial conditions lead to libration of Δω about 0° with a mean semiamplitude around 80°, whereas the rest of the sampled initial conditions seem to exhibit a mean semiamplitude around 180°, i.e., circulation (see Section 5.1 for more explanation of this dynamical behavior).

Zoom In Zoom Out Reset image size

Figure 7 shows an example of the 1000 yr extent of the orbital evolution simulation started from the best fit (i.e., maximum $-\mathrm{ln}{ \mathcal L }$; see Table 4). We show the evolution of the eccentricities, mutual period ratio P_rat, eccentricities e_b and e_c, and orbital inclinations i_b and i_c. The TOI-2525 system is consistent with moderate eccentricity evolution and appears to osculate outside of the low-order eccentricity-type 2:1 MMR. It is interesting that the secular evolution timescales are rather short, of the order of ∼120 yr. Therefore, future observations might be sensitive to transit depth variations. Furthermore, TOI-2525b may soon become a nontransiting planet (different timescales are possible within the mutually inclined posteriors).

Zoom In Zoom Out Reset image size

Our numerical orbital analysis of the system's configuration revealed that the TOI-2525 pair of planets is outside of the 2:1 MMR. However, the posterior of the apsidal alignment angle Δω shows bimodality with approximately equal fractions exhibiting libration and circulation. Thus, we took a closer look at the dynamical evolution of these two populations. The top and bottom panels of Figure 8 show the characteristic evolution of Δω for the same stable configuration as in Figure 7 (i.e., our best TTV+RV fit) and a random posterior fit with Δω in circulation, respectively. The left panels of Figure 8 show Δω as a function of time, whereas the right panels show the trajectory in the polar plot of ${e}_{b}{e}_{c}\cos {\rm{\Delta }}\omega$ versus ${e}_{b}{e}_{c}\sin {\rm{\Delta }}\omega$. In both cases, the polar plot shows that the system circulates around a point along the positive ${e}_{b}{e}_{c}\cos {\rm{\Delta }}\omega$ axis, i.e., an aligned configuration. For the best fit shown in the top panels, the trajectory in the polar plot is just small enough to miss the origin, and Δω exhibits large-amplitude libration about 0°. For the example shown in the bottom panels, the trajectory in the polar plot is just large enough to touch the origin, which leads to brief episodes of circulation of Δω. Thus, the bimodality of the apsidal alignment angle is simply due to slightly different sizes of the trajectories in the secular polar variables.

Zoom In Zoom Out Reset image size

We ruled out a 2:1 MMR librating configuration of the system based on the orbital posterior probability distribution constructed from the available transit and RV data. Similar to the TOI-2202 system, the osculating period ratio of the TOI-2525 pair of planets is well above 2, consistent with the prominent peak of period ratios of planet pairs observed by the Kepler mission (Lissauer et al. 2011; Fabrycky et al. 2014). Figure 9 shows an analytical analysis of the resonant and near-resonant dynamics in the 2:1 commensurability, following Nesvorný & Vokrouhlický (2016). The constant δ is an orbital invariant that defines the position of the system relative to the 2:1 period ratio, ψ is a combination of the resonant angles θ₁ and θ₂, and the variable Ψ is a combination of planetary masses, semimajor axes, and eccentricities (see Nesvorný & Vokrouhlický 2016 for details). Figure 9 is an updated version of Figure A1 of Trifonov et al. (2021), which now includes the position of TOI-2525, in addition to TOI-2202 (Trifonov et al. 2021), TOI-216 (Dawson et al. 2021), and Kepler-88 (Nesvorný et al. 2013). The median posterior probability values of TOI-2525, listed in Table 4, lead to δ ≃ −1.60; therefore, the system is firmly outside the libration region, together with the Kepler-88 and TOI-2525 systems. The only system that is in a 2:1 MMR is TOI-216 (see Dawson et al. 2021 for details).

Zoom In Zoom Out Reset image size

An interesting question is how frequently state-of-the-art planet formation models produce systems like TOI-2525, TOI-2202, or TOI-216 and what drives their formation. We explored the abundance and origins of systems with multiple giant planets in synthetic planet populations from the Generation III Bern model (Emsenhuber et al. 2021a; Schlecker et al. 2021a; Emsenhuber et al. 2021b; Schlecker et al. 2021b; Burn et al. 2021; Mishra et al. 2021). The model includes the mechanisms relevant for the dynamical evolution of multiplanet systems, in particular types I and II planet migration (Paardekooper et al. 2011; Dittkrist et al. 2014), eccentricity and inclination damping through planet–disk interaction (Cresswell & Nelson 2008), and dynamical evolution modeled via an N-body integrator (Chambers 1999). The planet population with a 0.7 M_⊙ host star introduced in Burn et al. (2021) is most suitable for the comparison with the K dwarf system presented here. Out of its 999 systems, 92 contain planets more massive than 30 M_⊕. We find 51 systems with more than one such planet within orbital periods of 1000 days in the population, typically around host stars with enhanced metallicity. Only a single system includes a pair of warm giant planets in an MMR-like configuration; it contains three giant planets with periods of 48.0, 96.5, and 479.8 days.

The system emerged from a numerical disk with a large solid material content, which is reflected in a high dust-to-gas ratio of 0.024 as compared to the population median of 0.014. This led to efficient core growth and runaway gas accretion of three protoplanets. Through simultaneous inward migration and gravitational interaction, the two warm gas giants were eventually captured in an MMR, a configuration that persisted until the end of the N-body integration at 20 Myr. This comparison proposes that the TOI-2525 system was formed in a similar metal-enriched disk. Thus, TOI-2525 is in a configuration whose realization through core accretion is rare but possible, according to the simulations.

We note that the literature contains other examples of pairs of warm Jovian planets that have been extensively studied, for instance, the Kepler-9 system, consisting of a G dwarf star orbited by a mini-Saturn planet pair in a 2:1 MMR (Holman et al. 2010). Other examples include the HD 82943 (Tan et al. 2013) and TIC 279401253 (Bozhilov et al. 2023) systems, both of which are G dwarf stars with almost identical 2:1 MMR Jovian-mass planet pairs. Although rare, M dwarfs have also been found to host 2:1 MMR warm massive systems, with the GJ 876 multiple-planet system (Rivera et al. 2010; Millholland et al. 2018; Trifonov et al. 2018) being widely considered as a benchmark for planetary dynamics and planet formation. The formation of such systems may not depend on the stellar type, but its occurrence rate remains an important observable to be further studied.

Objects TOI-2525b and c are two warm, low-density giant planets, a category of planetary systems whose frequency is steadily growing in the literature. The estimated mean density of TOI-2525c of ρ_c = ${1.014}_{-0.076}^{+0.084}$ g cm⁻³ is lower than that of Jupiter (ρ_Jup = 1.33 g cm⁻³) but higher than that of Saturn (ρ_Sat = 0.69 g cm⁻³). Therefore, the density of the Jovian-mass planet TOI-2525c is not surprising. The Neptune-mass planet TOI-2525b, however, has a mean density of ρ_b = ${0.174}_{-0.015}^{+0.016}$ g cm⁻³, which makes it among the lowest-density Neptune-mass planets known to date. The density of TOI-2525b is a bit larger than the Kepler-51 b, c, and d "super puffs" (Steffen et al. 2013; Masuda 2014) and is comparable to low-density planets like WASP-107 b (Rubenzahl et al. 2021), WASP-131 b, WASP-139 b (Hellier et al. 2017), WASP-21 b (Barros et al. 2011; Ciceri et al. 2013), and HATS-46b b (Brahm et al. 2018), among others. We use the mass–radius models of Fortney et al. (2007) to infer the core mass of planets TOI-2525b and c given the estimated masses, radii, and stellar parameters. We calculate core masses of 6.9 ± 1.4 and 36.6 ± 9.3 M_⊕ for TOI-2525b and c, respectively. The ratio of core mass to total mass Zp is therefore Zp_b = 0.25 ± 0.05 and Zp_c = 0.17 ± 0.04.

Figure 10 shows a mass–radius plot for all known transiting planets with measured masses validated by TTVs or RVs. Figure 10 is limited to the giant planets in the range of 0.05–11 M_Jup and 0.15–3.0 R_Jup. Panel (a) of Figure 10 is color-coded with the planetary equilibrium temperature (T_eq), whereas panel (b) is color-coded with the planetary density. In both panels, we plot an interpolated model for the mass–radius relationship assuming the estimated luminosity of TOI-2525, semimajor axis = 0.2 au, and an age of 3.1 Gyr from Fortney et al. (2007). From Figure 10, it is clear that the higher the T_eq of the planet, the larger the radius. Except for the massive planets with a few Jupiter masses or more, the larger radius correlates with low planetary density. Object TOI-2525c is consistent with the mass–radius model from Fortney et al. (2007); TOI-2525b is consistent with a large radius, and its mean density is among the lowest for the Neptune-mass planets.

Zoom In Zoom Out Reset image size

We also model the evolution of both planets in the system, using CEPAM (Guillot & Morel 1995; Guillot et al. 2006) and a nongray atmosphere (Parmentier et al. 2015), to provide constraints on their interiors. We assume simple structures consisting of a central dense core surrounded by a hydrogen and helium envelope of solar composition. The core is assumed to be made of 50% ices and 50% rocks. Figure 11 shows the resulting evolution models. The core mass of TOI-2525c is found to be between 20 and 43 M_⊕. This indicates that the enrichment in heavy elements of TOI-2525c could be comparable to Jupiter's, which is between 8 and 46 M_⊕ (Guillot et al. 2022). With a radius similar to Saturn's (1.08 times larger) but 3.4 times less massive than Saturn, TOI-2525b is an uncommon example of a very low density and inflated planet with an equilibrium temperature close to 500 K. The H–He envelope of TOI-2525b is found to be between 19 and 24 M_⊕. The case of TOI-2525b is challenging for the traditional core accretion formation scenario. With our simple modeling of TOI-2525b, such a small envelope hints that the accretion of H–He has potentially been hindered. Characterizing the atmospheres of both planets of the system would be very useful to understand their structure and formation.

Zoom In Zoom Out Reset image size

Not many inflated Neptune-mass planets are known, making TOI-2525b a useful addition to the sample of transiting planets with a measured mass. The low density, and thus large scale height, of TOI-2525b makes it a good target for a future atmospheric investigation with transmission spectroscopy.